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1 Introduction 



In the companion paper [T], we proposed a generalization of the Petviashvih iteration 
method for finding stationary soUtary waves m(x) of scalar and vector Hamiltonian equations 
with arbitrary form of nonlinearity: 



where M is a self-adjoint differential operator and, in the vector case, the nonlinear term 
must satisfy a condition dFi/duj = dFj/dui. (Recall that the original Petviashvih method 
[2] was proposed for scalar equations with power-law nonlinearity F{x, u) = vP.) A common 
form of operator M (in the scalar case) is 



where /x is the propagation constant of the solitary wave. Thus, the generalized Petviashvih 
method, that obtains solutions with a specified propagation constant, can be applied to the 
same class of equations as the well-known imaginary-time evolution method (ITEM) (see, 
e.g., [3]-[6]) that is used to find solitary waves with a specified power. 

In the present work we extend the results of [I] as follows. In Section 2, we establish 
a mathematical relation between the generalized Petviashvih method and the ITEM. This 
discussion will also set the stage for the main result of this work, presented in Section 3. 
There, we develop the ideas behind the original and generalized Petviashvih methods [71 [T] 
and propose a new technique that we refer to as the mode elimination. This technique 
can be used to obtain nonfundamental (see below) solitary waves, which the methods of 
[l]-[^ cannot obtain (the iterations would diverge). However, since alternative methods 
of obtaining nonfundamental solitary waves exist [3l [8] , we see the main use of the mode 
elimination in that it can considerably accelerate convergence of various iteration methods. 
The corresponding examples are presented in Section 4, and the summary of our results is 
given in Section 5. 

2 Convergence rates of the Petviashvili and the imaginary- 
time evolution methods 

In this Section, we will compare the convergence properties of the generalized Petviashvili 
method [l] with those of the accelerated ITEM proposed in Ref. [6]. This discussion will 
highlight a feature of the generalized Petviashvili iteration scheme that will be important 
when we present our main result — the mode elimination technique — in the next Section. 
Of the two versions of the ITEM (with power and amplitude normalizations) considered 
in [6], we will focus on the one with power normalization, because its linearized operator 
can be readily compared with that of the generalized Petviashvili method. In order not to 






M = /i- V^ 
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obscure the main ideas by technical details, we restrict our presentation to the case of a 
single real- valued equation (II. Ih with M given by Eq. (|1.2p . i.e., to: 

V2n + F(x,u) = ^u. (2.1) 

It is well-known that the convergence of an iteration method is determined by the 
properties of the linearized iteration equation. Namely, let n„ be the solution obtained at 
the nth iteration, and let the "error" Un be defined as 

Un = Un — \Uri\ ^ \u\ . (2-2) 

As will be shown below, it satisfies a linearized iteration equation of the form 

= (1 + Ar/:)n„, Ar>0, (2.3) 

where C is the linear operator that results when the iteration method is linearized on the 
background of the solitary wave u, and At is an auxiliary scaling parameter. From a 
conceptual point of view, the presence of Ar emphasizes the analogy of iteration methods 
with numerical methods of solving time-dependent differential equations (see, e.g., 19] ); 
from a practical point of view, it can be used to ensure (in certain cases) or optimize the 
convergence of the method, as we will discuss later on. 

Let us begin with general remarks regarding the convergence rate of the linearized 
iteration equation (|2.3p . Suppose that the eigenfunctions of C form a complete set in an 
appropriate functional space, so that Un can be expanded over them. Let the minimum and 
maximum eigenvalues of C be Amin and Amax- Then the convergence rate of the iteration 
method can be defined as log(l/i?), where the convergence factor R is the maximum (in 
magnitude) eigenvalue of the operator on the r.h.s. of (|2.3p : 

i? = max{|l + AmaxATl, |l + AminAr|}. (2.4) 

Clearly, R <1 needs to hold in order for the iterations to converge, which implies 

Amax < and 1 + A^inAr > -1 . (2.5) 

Moreover, if Amax = 0, then the corresponding eigenfunction of C needs to be a translational 
eigenmode (if it exists) of the linearized Eq. (12. ip . which only shifts the solution in space 
and hence does not affect the convergence of the method. The smaller the convergence 
factor R, the faster the convergence. It can be readily shown [6] that the minimum value 
of R occurs at 

An = ^ . (2.6) 

■"-mm \ -"-max 

(recall that Amin < Amax < 0) and equals 

_ 1 ~ (Amax /Amin) 
1 + (Amax/A mm^ 
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Therefore, the closer the ratio (Amax/Amin) to 1, the faster the convergence of the itera- 
tion method. Below we will compare the possible values of this ratio for the generalized 
Petviashvili method and the accelerated ITEM. To that end, we first need to cast the 
linearizations of these methods into the form of Eq. (j2.3p . 

Let Lq denote the nonlinear operator of Eq. (j2.1|) (or, more generally, of the stationary 
wave equation whose solution we are looking for), so that that equation is rewritten as 

Lon = 0. (2.8) 

Let L be the corresponding linearized operator, so that 

Lq{u + u) = Lqu + Lu = Lu , for any \u\ <^ \u\ . (2-9) 

Note that for Hamiltonian wave equations, L is always self-adjoint. With these notations, 
the generalized Petviashvili method is [1]: 

Un+l-Un=\N {LoU)n-J—, Un ] At , (2.10) 

V {Un,NUn) J 

where 

7 = 1 + ^. (2.11) 

Here and below, the inner product between two real- valued functions is defined in a standard 
way: 

/oo 
/(x)5(x)dx. 

For the positive definite and self-adjoint operator N in (j2.10p . we take the simplest form 
used in [1]: 

N = c-V^, (2.12) 

where the constant c is given by Eq. (3.11) of [Ij. The constant a in Eq. (j2.1ip above is 
such that 

Lu aNu (2.13) 

in a certain least-square sense; a formula for computing this constant at each iteration can 
be found in either of Eqs. (3.12) or (3.15) of [Ij, but will not be needed here. Following 
the steps of a calculation found at the beginning of Section 2 of [1], it is straightforward to 
show that the linearized form of the generalized Petviashvili method is: 

un+i -un= [N-^Lun - 7^^^ ^ ) Ar . (2.14) 



Next, the accelerated ITEM of Ref. [6] is: 



Un+l 



P 



Un+l, (2.15) 



XUn+l,Un+l) . 

Un+l -Un = ( V^n„ + F(x, Un) - fJ-nUnj Ar, (2-16) 
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Un) 



(2.17) 



where P = u^d^ is the specified power of the solitary wave. The positive definite and 
self-adjoint operator K is referred to as the acceleration operator for the ITEM [3l[6]. For 
simplicity, we take K to have the same form (j2.12p as the operator N in the generalized 
Petviashvili method, with the c being now an arbitrary positive constant. The linearized 



Thus, the "primordial" operator in the linearized equations of both the generalized Petvi- 
ashvili method and the accelerated ITEM has the form: 



With L being the linearized operator of (j2.ip . the continuous spectrum of L is an interval 
(or, when F{x, u) is a periodic function of x, a union of intervals), one of the end points of 
which is A = —1 (see, e.g., |6J and references therein). This eigenvalue of L corresponds to 
the eigenvalue A = — oo of L. Then a possible spectrum of L is shown in Fig. [1^. 

Even though the first terms on the r.h.s.'s of ()2.14p and ()2.18p have the same form 
(j2.19p . the eigenvalues of the corresponding operators C are different for two reasons. First, 
the values of c in operators and K are, in general, different, which makes different the 
eigenvalues of the corresponding L's. Second, the nonlocal terms (involving inner products) 
in (|2.14p and (j2.18p modify the eigenvalues of L in different ways. We now consider this 
latter issue in more detail. 

In regards to the operator of the linearized Petviashvili method (j2.14p . we recall a fact 
[Ij that is important for our discussion both here and in the next Section. Namely, the role 
of the nonlocal term in that operator is to (nearly) eliminate from Un+i the eigenfunction 
of L = N~^L whose profile is close to that of the solitary wave u, while leaving the other 
eigenfunctions and their eigenvalues (nearly) unchanged. This is ensured by taking the 
constant 7 and the operator N to satisfy (12. lip and ()2.12p , respectively. The adverb "nearly" 
is used above to account for the fact that relation (12.13P for Eq. (j2.ip with a general 
nonlinear function F{x, u) holds only approximately. It is exact only for wave equations 
with power-law nonlinearity [7J, for which the original Petviashvili method was proposed [2]. 
However, the special choice of the constant c in (12.12p . as noted after that equation, makes 
the approximation in (j2.13p sufficiently accurate at least near the "core" of the solitary 
wave. 

Continuing with the discussion about the effect of the nonlocal term in p.l4p on the 
eigenvalues of the corresponding operator C, let us suppose that ti is a fundamental solution 
of the nonlinear wave equation. (E.g., in the case of Eq. ()2.ip . the fundamental solution. 



form of ITEM (|2T5D - (|2T7]) is [6J: 




(2.18) 



L = (c- V2)-1l. 



(2.19) 
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unlike nonfundamental ones, has no nodeqj. For a more general Eq. (jl.ip where the 
operator M is different from V^, fundamental solutions may have nodes (as, e.g., the lump 
solution of the Kadomtsev-Petviashvili equation [10]); in that case, their distinguishing 
feature is that they have one "main" hump, while the nonfundamental solutions usually 
have several "main" humps.) Then the "u-like" eigenfunction of operator N~^L mentioned 
in the previous paragraph (see also ()2.13p ) corresponds to the largest eigenvalue, Ai, of that 
operator; see Fig. [1^. Since this eigenfunction is eliminated by the nonlocal term at each 
iteration, then the resulting spectrum of the operator on the r.h.s. of (|2.14p is as shown 
in Fig. [Da. Thus, for this operator, Amax ~ ^2 and Amin ~ Amin; the reason for using 
"~" instead of "=" is that relation (j2.13p holds approximately, as we noted above. Now, if 
A2 < and the step size At satisfies a condition 



then according to (j2.5p . the generalized Petviashvili method converges to u. As a sidenote, 
we mention that for equations with power-law nonlinearity, L is known [llj to have only 
one positive eigenvalue, and hence the Sylvester inertia law (see, e.g., Theorem 7.6.3 in [12j) 
guarantees that Ai is the only positive eigenvalue of L = N~^L. 

Now let us consider the hnearized operator C in (f238]) for the ITEM (f2T5]) - (f2T7ll . In 
[6], we showed that the set of discrete eigenvalues of this C is the union of two sets: (i) the 
roots of a function 



where 'i/'j is the eigenfunction of L corresponding to the eigenvalue Xj, and also (ii) the 
set of those Xj for which (u, ipj) = 0. This is shown schematically in Fig. [it, with ■03 there 
satisfying (u, ip^) = 0. (Note that Q(A) does not need to be defined for the continuum 
eigenvalues A.) Thus, for the operator C in (j2.18p . Amin > Amin and Amax > A2. 

The consideration of the two preceding paragraphs shows that even when the acceler- 
ation operators N and K in (j2.14p and (j2.18p are the same (i.e., have the same c), one 
cannot, in general, make a definite statement on whether the ratio (Amax/Amin)! and hence 
the convergence rate, is greater for the generalized Petviashvili method or for the acceler- 
ated ITEM. Moreover, the fact that the values of c in and K are generally different, 
and hence so are the eigenvalues Xj of the corresponding two L's, further obstructs the 
comparison of the convergence rates of the two methods. The only two statements that can 
be made here are the following, (i) For equations (j2.ip with arbitrary nonlinearity, if the 
ITEM converges to a fundamental solution, then we expect that in most cases (see below), 
so does the generalized Petviashvili method, (ii) For equations (j2.ip with power-law non- 
linearity -F(x, u) = vP, the Petviashvili method with the optimal choice of Ar converges to 
the fundamental solution faster than does the optimally accelerated ITEM p.l5p - ()2.17p . 
^By nodes in D > 1 spatial dimensions, we mean sets of points of dimension less than D where u(x) — 0. 



l-FAminAr > -1 



(2.20) 




(2.21) 
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To justify statement (i), first recall that for fundamental solitary waves, 

(Amax)petviashvili ~ A2 < (Amax)lTEM j (2.22) 

as long as the value c in the operator (j2.19p is taken to be the same for both methods. 
Next, if the ITEM converges, then according to (j2.5p . (Ajnax)iTEM < 0, thereby implying 
that A2 < 0. However, by the Sylvester inertia law, the sign of A2 does not depend on 
the actual value of c (as long as c > 0). Therefore, with a possible exception of those 
cases where A2 is close to zero, the left part of (|2.22p yields (Amax)petviashviii < 0, which 
means that the generalized Petviashvili method converges. To prove statement (ii), first 
note that operator L in this case satisfies the conditions of Theorem 4 of Ref. [6j, so that 
c = ;U is the optimal value for K and (Amin)iTEM = Amin (= — !)• Next, in the Petviashvili 
method for the equation with F{x, u) = u^, = M [21 |T] and hence c = as well, whence 
(Amin) Petviashvili = Amin- Thus, in this casc, 

(Amin) Petviashvili — (Amin )lTEM- (2.23) 

Combining Eq. (j2.23p and inequality ()2.22|) . where now the sign "~" must be replaced with 
"=", one concludes that (Amax/Amin) should be greater for the Petviashvili method; hence 
statement (ii) follows. 

A simple example illustrating statement (ii) is the stationary nonlinear Schrodinger 
equation in one dimension: 

Uxx + = u, \u\ as |x| 00, (2.24) 

for which the ITEM with the parameters c = //(= 1) and At = 1.5, corresponding to the 
optimal acceleration, converges to the accuracy of 10~^'^ in 33 iterations. The Petviashvili 
method (|2.10p with At = 1.5, alpha = 2 (as in the original Petviashvili method; see [1]), 

and 7 given by ()2.1ip . converges to the same accuracy in 19 iterations. Here both methods 

2 

start with the initial condition uq = . In our numerical experiments of finding the 
fundamental solutions of non-power-law equations (not covered by the above statement 
(ii)), we also observed that the generalized Petviashvili method is faster than the optimally 
accelerated ITEM (l2T5]l - (f2Tni : see, e.g.. Example 3.1 in [I]. (The ITEM with amplitude 
normalization [6] can still be faster than the generalized Petviashvili method.) 

However, in a situation where both methods converge to a nonfundamental solitary wave, 
the optimally accelerated ITEM can be faster than the generalized Petviashvili method. As 
an example, let us revisit the equation with a double-well potential: 

Uxx + V{x)u-u^ = ^u, y(x) = 6 (sech^(x - 1) + sech^(x -h 1)) , (2.25) 

considered in Example 3.2 of [T]. We will focus on its anti-symmetric solution (see Fig. [2^) 
with the propagation constant /i = 1.43 and the corresponding power P = dx = 10. 
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This solution is nonfundamental since it has a node; the fundamental solution in this case 
is a two-humped pulse with its maxima located near the maxima of the potential. The solid 
and dashed lines in Fig. [2]3 show the evolutions of the error norm, defined as 

V {Un,Un) J 

for the generalized Petviashvili method and the optimally accelerated ITEM, respectively. 
In both cases, the parameter At was emprically optimized (see (j2.6p ) to yield the maximum 
convergence rates; the respective values are Ar^, ^Petviashvili = 1-6 and At^kjtem = 0.7. Also, 
in the case of the generalized Petviashvili method, the value c = 5.04 was algorithmically 
computed [Ij, while for the ITEM, c = 1.5 was empirically found to yield the optimal 

2 

convergence rate. As the initial condition for both these methods, we took uq = 2xe~^ . 
As seen from Fig. [2]3, the optimally accelerated ITEM is about one and a half times faster 
than the generalized Petviashvili method. The reason behind this can be understood by 
looking at the spectra of the corresponding operators L in (j2.19p with the above values 
c = 1.5 for the accelerated ITEM (Fig. [3h)) and c = 5.04 for the generalized Petviashvili 
method (Fig. [3t). Namely, when one starts with an anti-symmetric initial condition (as 
we did above), the symmetric eigenmodes corresponding to X2k+ii k = 0,1,... do not 
contribute to the error tt„. Then from (j2.7p and Figs. [3)3, c, 

jj ^ 1 - (A2, c=i.5/Amin, 0=1.5) 1 - 0-41 . 

1 — 1-^2, c=1.5/Amin, c=1.5j i + U.4i 

jj ^ (Amax continuum, c=5/Amin, c=5) ^ 0.28 

-"-Petviashvili ~ -] i ( \ 7\ \ ~ i , n 98 ~ "J-OD) 

~r l,^max continuum, c=5/ ^min, c=5 J ^ i o.^o 

and hence the corresponding numbers of iterations to reach the accuracy of 10^^*^ can be 
estimated as: 

-10 In 10 -10 In 10 

'T-max, ITEM ~ ^ — 26, nmax, Petviashvili ~ 1 7^ — 40. 

In H ITEM hi H Petviashvili 

These estimates are in very good agreement with the numbers of iterations (25 and 37, 
respectively) reported in Fig. [2)3. Note also that the empirically found optimal values of 
At* reported above agree with Eq. ()2.6p and the spectra shown in Figs. [3b,c. 

3 Mode elimination technique for improving convergence of 
iteration methods 

Here we develop the ideas of Ref. [T] and extend the generalized Petviashvili method so 
that it could be employed for two additional purposes: (i) obtaining certain nonfunda- 
mental solutions of stationary nonlinear wave equations; and (ii) accelerating convergence 
of iterations methods. We emphasize that the technique we propose can be applied to 
any iteration method and to single and coupled equations as well. For simplicity of the 
presentation, below we illustrate it for single equations of the form (j2.ip . 
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We begin with the observation that in most cases (with Eq. (j2.25p being a notable 
exception), the generahzed Petviashvih method would not converge to a nonfundamental 
solution of a given wave equation. The reason for that can be understood from the following 
simple example. Consider an equation 

Uxx + (6sech^x + u'^)u = fiu . (3.1) 

When the amplitude of u is small, (13. ip has two solutions: the fundamental, {u^^^ ~ 
esech^x, /i*-^-* ~ 4}, and the nonfundamental, {u^"^^ ~ esechxtanhx, /i*-^-* ~ 1}) where 
e <C 1. Then the operator obtained by the linearization of Eq. (j3.ip on the background of 
the nonfundamental solution, 

L^dl + Gsech^x - /i^^^ , (3.2) 

has two largest eigenvalues: Ai k, ijS^^ — ji^"^^ ss 3 > and A2 ~ n^'^^ — fJ^^'^^ = 0, with the cor- 
responding eigenfunctions being approximately u^^^ and u^'^\ As we noted in Section 2, the 
nonlocal term in the linearized iteration equation (j2.14p nearly eliminates the eigenfunction 
of operator L = N~^L which is "similar" to the background solution u^'^h However, the 
eigenfunction of L corresponding to the eigenvalue Ai > of L is not eliminated, and hence, 
according to the discussion found before Eq. (j2.20p . the generalized Petviashvili method 
will not converge to solution u^'^\ 

The above example suggests a simple way in which the generalized Petviashvili method 
(j2.10p can be modified so that it would converge to a nonfundamental solution u (given, of 
course, an initial condition close to u). In the general form, this modified method is 



Ar, 



Un+l U-n — JV \^Q"')n 1 , at \ 'unst , fj) , 'Puns' 

(3.3) 

where 7 and N are defined as in ()2.1ip and (I2.12|) . ^^nst ^^'^ the functions that approximate 
the eigenmodes of operator (N^^L) with positive eigenvalues (excluding the background 
solution u), Junst is the number of such eigenmodes, and 

/unst -"^ ' (i\ . ! "unst , a) (,;) , • \"-^y 

Here a^^^^t^ defined analogously to (j2.13p : 

L<p\iL « (^^£tN<p\iL , (3.5) 

is computed according to Eq. (3.12) of [1]. In the context of the example in the previous 
paragraph, J^nst = 1 and (j)^^^^^ = u^^^ . 

Following the lines of the analysis of Section 2 in Ref . [1] , it is straightforward to show 
that in method (|3.3p , ()3.4p , the components of the error "aligned along" the modes (pl^^^i^ 
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j = 1, . . . , Junsti are nearly eliminated at every iteration; this is guaranteed by the form of 
the coefficients 7unsf Therefore, in what follows, we refer to method (j3.3p as the mode 
elimination method. In Section 4 below, we will present the results of applying this method 
to a two-dimensional equation of the form (j3.ip to obtain its nonfundamental solutions. 
Remark It is clear that the success of the mode elimination method hinges upon the 

knowledge of the "unstable" eigenmodes (p^^^st- However, in many cases, an approximate 

(7) 

knowledge of (/>unst suffice. 

We now show how the mode elimination technique can be used to accelerate convergence 
of iteration methods. The reason that a given method converges slowly is, according to (j2.7p . 
that the ratio Amax/Amin is small. Since for an appropriately chosen operator N, |Amin| = 
0(1) (see, e.g.. Figs. [3l3,c), then for a slowly convergent method, the eigenvalue lAmax must 
be small. Then if one can eliminate the corresponding eigenmode, similarly to how it is 
done in (jS.Sp . one essentially replaces (Aniax)oid with (Ajnax)new < (Amax)oid 

(< 0). Then 

the ratio Ajnax/Amin increases and so does the convergence rate of the iteration method. 
The practical issue here is how to find the mode, 0siow) which slows down the convergence. 
Fortunately, this is rather easy to do using the following observation. For At < Ar*, where 
At* is defined in p.6p . the factor (1 + AsiowAr) = (1 + AmaxAr), which governs the decay 
of 6 

siow! is the largest among such factors for all the eigenmodes of [N ^L). Then after 
some iterations, the content of the error Un = Un — u becomes dominated by the eigenmode 
0siow) and hence 

(/•slow oc (u„ - ■u„_i) . (3.6) 

The elimination of the function (u„ — Un-i) is carried out in exactly the same way as in 
(j3.3p . yielding the method: 



Un+l 



where 



,r-Ur N {Lou)n) (0slow,n, (i'0'")n) , 
iV [LoU)n - 7 —/ ^7 ^ Un - 7slow, n JJ TTJ T 0slow, n 



Ar, 
(3.7) 



V'slow, 71 — Un—li 7slow,7i — J- H T , Onslow, 71 — TT 777 T ' v"^'"/ 

Onslow, 71^''" \'Pslow,n) -' ' 'Pslow, n/ 

Note the coefficient s in (j3.8p . which we will comment on in the next paragraph. We will 
also provide examples that demonstrate the efficiency of the accelerated Petviashvili method 
(j3.7p . ()3.8p and its extensions to other iteration methods, in the next Section. 

Similarly to the analysis of Ref . [1] , one can show that the role of coefficient s in (j3.8p 
is to control how much of the mode (/>siow,n is subtracted at each iteration. We found 
empirically that in most cases, it is beneficial for the convergence rate to subtract not the 
entire i;^siow,n-component from Un but only part of it, usually somewhere between 40% and 
80% (i.e., use s ~ 0.4-0.8). (However, even using the value s = 1 leads to a significant 
increase in convergence rate compared to the corresponding non-accelerated method when 
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the latter is slow.) The justification of using < s < 1 (or, alternatively, 1 < s < 2) rather 
than s = 1, is based on the same considerations, found before Eq. (13. 6p . which led us to 
propose the accelerated method ()3.7p . Namely, to uphold those considerations, (j)siow,n is 
to remain the most slowly decaying eigenmode of (N^^L) at every iteration. In the case 
where the entire amount of it is subtracted at the (n + l)st iteration, it is not obvious (and 
probably not true) that the error Un+2 at the next iteration would consist mainly of the 
mode (/>siow, n+i = Un+i — Un, which will be subtracted at the (n + 2)nd iteration. However, 
if only s ■ 100% of mode c/'siow, n is subtracted, this mode can still remain the most slowly 
decaying as long as 

1(1 - s) • (1 + AsiowAr)! > |1 + AnextArl, (3.9) 

where Asiow is the eigenvalue corresponding to (j)siow,n, and Anext is the eigenvalue corre- 
sponding to the next most slowly decaying mode. Yet, for s not too small, the l.h.s. of (13. 9p 
is still considerably less than |1 + AsiowAr|, and hence the convergence rate of the original 
iteration method is increased. 

To conclude this Section, we compare our mode elimination technique for convergence 
acceleration with the Steffensen's method (see, e.g., [13]), which is based on applying the 
Aitken's acceleration algorithm every given number of iterations. The idea of the Stef- 
fensen's method is the following. Suppose one has three consecutive iterative solutions Un, 
Un+i, Un+2 about which one knows that they satisfy 

« for all X, (3.10) 

u„+i(x) u,„(x) 

where Un is the error defined in l\2.2\) . Using these solutions, one applies the Aitken's 
algorithm: 

_ A {Un+l — Un)"^ /o n\ 

Un+3 = U^=Un ■ , (3.11) 

Un+2 - ^Un+l + Un 

and then proceeds to computing the next few iterations Un+A, ... ,Un+n^^cci+2 with the 
original iteration method, where naccei > 3. Then one uses Un+n^^^^i, Un+n^^^^i+i, Un+n^^^^i+2 
to compute "Un+naccoi (|3.1ip with n — > n + riaccel, and so on. In [H], this method was 
successfully used to accelerate the convergence of the original Petviashvili method for the 
nonlinear Schrodinger equation in 3 spatial dimensions. 

Aitken's algorithm (j3.1ip systematically reduces the error {u^ — u) only when (j3.10p 
holds sufficiently well, which occurs under the same condition (j3.6p that must hold in order 
for the mode elimination method to work. However, the sense in which (j3.6p is to hold 
is drastically different for these two acceleration techniques. For the mode elimination, it 
suffices if (j3.6p holds approximately near the "core" of the solitary wave, since i;^siow enters 
Eqs. (j3.7p . (j3.8p via the inner products with functions that are essentially nonzero only 
in that spatial region. On the contrary, for the Steffensen's method, p.lOp has to hold 
pointwise and, in particular, far away from the "core" of n(x). In the latter spatial region, 
the denominator of ()3.1ip is nearly zero, and hence even a small ripple in n„, Un+i, or 
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Un+2 can result in a large distortion of u„ . This was indeed observed in our numerical 
experiments, except in the cases where |Amax| *C |Ancxt|, where |Ancxt| is defined after 
()3.9p . Thus, we expect our mode elimination technique and the Steffensen's method to 
be competitive in those latter cases, but expect the mode elimination technique to have 
superior performance over that of the Steffensen's method when there are more than one 
eigenmodes with A ?s Amax- This expectation is borne out by Examples 4.2 and 4.3 reported 
below. 



4 Examples of the mode elimination technique 

Below we illustrate the application of the mode elimination technique to obtaining nonfun- 
damental solitary waves and to accelerating convergence of iteration methods for Eq. ()2.ip . 
In Ref. [8], we already showed by extensive simulations that this technique can greatly 
accelerate convergence of a class of universally-convergent iteration methods for both sin- 
gle and coupled equations. (Method (j4.4p presented below is a particular member of that 
class.) Therefore, here we will focus on clarifying the role of parameter s in Eq. (13.80 for 
optimizing the convergence rate and also on demonstrating the applicability of the mode 
elimination technique to various classes of iteration methods. 

Example 4.1 Here we will demonstrate that method (|3.3p . (j3.4p can be used to obtain 
nonfundamental solitary waves when approximate information about the unstable eigen- 
modes of (N^^L) is available. We will also compare the performance of this method with 
that of a universally-convergent method proposed in [8] . 

Equation 

V'^u + Vo{sech X sech y)'^ u + = IJ.U, Vq = 20 (4.1) 



is a two-dimensional counterpart of Eq. ()3.ip . Since the potential well in (14. ip is sufficiently 
deep (Vq ^ 1), this equation admits several nonfundamental solutions. Below we report the 
details of finding the first of them which corresponds to /i = 8 and is shown in Fig. [H For 
this solution, we expect the generalized Petviashvili method to have one unstable eigenmode 
(in addition to the mode approximated by u that may possibly also be unstable), and 
approximate this eigenmode by 



•'unst 



e-^(^/^)', r^ = x^ + y\ (4.2) 



The width W in ()4.2p is found iteratively from the formula 

^2^2(n^^^^^ (4.3) 

3 {Un, Un) 

in deriving which we assumed that u (x x (punst ■ Starting with the initial condition uq = 
2x e~^^^^^^\ method (j3.3p . (j3.4p with a nearly optimal At = 0.7 took about 50 iterations 
to reach the accuracy of lO^^*'. Thus, the generalized Petviashvili method with mode 
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elimination (j3.3p . (j3.4p converges to this nonfundamental solution, while the generalized 
Petviashvili method ()2.10p without the mode elimination diverges. 

We also obtained the same solution by a method based on the "squared" operator 

Un+1-Un = - [N LN LoU)n-Tn -, ^7 ^ Ar. (4.4) 

[ (m„, NUn) J 

(The name "squared" comes from the fact that (N-'^L)'^ appears in the linearized version 
of (j4.4|) .) In [8], we showed that this method belongs to a family of universally-convergent 
methods (i.e., methods which can converge to any nonfundamental solution of a given 
equation provided that the initial condition is sufficiently close to that solution) for either 
of the following choices of r„: r„ = or 

" ^ " (K, iLN-^Lu)n) / {un, Nur.)) At ■ ^^-^^ 
Note that this r„ is defined similarly to 7„ in the generalized Petviashvili method (see Eqs. 
()2.1ip and (j3.4p ). Since we are looking for a nonfundamental solution of (j4.ip . then using 
the value for r„ given by (j4.5p as opposed to r„ = will not eliminate the mode with 
the maximum eigenvalue (see the discussion after Eq. (j3.2p ). and hence will not speed up 
the convergence of the iterations. Therefore, in the remainder of this Example we report 
the results for method (14. 4p with = 0. Starting with the same initial condition as 
above, this method with the operator N computed as in [ij and with a nearly optimal 
Ar = 0.5 took about 190 iterations to converge to the accuracy of 10^^'^. Thus, the mode 
elimination method (j3.3p . (j3.4p is several times faster than the squared-operator method 
(|4.4p for finding the first nonfundamental solution of (j4.ip . (We also observed that method 
()3.3p . (|3.4p is less sensitive to the choice of initial conditions than method (j4.4p .) However, 
when we additionally included the step of eliminating the slow mode, as in Eqs. (j3.6p - 
(I3.8p . into both methods, the difference in their convergence rates was significantly reduced. 
Namely, the convergence of method (13. 3p . (|3.4p . which has already been quite rapid, was 
not improved by this additional step (and the number of iterations remained around 50), 
while the squared-operator method now took about 70 iterations to converge. 

We also applied both methods to finding the second nonfundamental solution of (j4.ip . 
which has the shape similar to ^4(1 — Br^) e~^^^^^^ with = + and A,B,C = const 
(see Fig. [7] below). For this solution we found, through experimentation, that one needs 
to include five unstable modes into (j3.3p . For the respective optimal Ar's, the generalized 
Petviashvili method with mode elimination (j3.3p . (13.41) was found to be about 50% faster 
than the squared-operator method (|4.4p . However, this advantage in the convergence rate 
is offset by the increased complexity arising from the need to guess the number and profiles 
of unstable modes and then to estimate their parameters (namely, the widths). Therefore, 
we conclude that the mode elimination method may be more efficient than the squared- 
operator method for finding the lowest- order nonfundamental solitary waves, as long as 
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some reasonable guess about the unstable modes can be made. However, for finding second- 
and higher-order nonfundamental solutions, method (j4.4p appears to be easier to implement 
and hence more practical. 

Example 4.2 In this and the next two Examples, we demonstrate the efficiency of the 
convergence acceleration technique based on the mode elimination, as in (j3.7p and (j3.8p . 
for three different iteration methods. In this Example, we apply this technique to the 
generalized Petviashvili method. 

We look for the fundamental solitary wave of an equation arising in the theory of non- 
linear photonic lattices: 

V^M -|- Vo(cos^ X + cos^ y)u + u'^ = fiu. (4-6) 

for three choices of the potential amplitude and the propagation constant: 

(a) : ^0 = 4, /i = 4.95; (6) : Vq = 4, ^l = 6.5; (c) : Vq = 0, /i = 1 . 

In case (a), the propagation constant is close to the edge of the continuous spectrum band, 
and the solitary wave occipies many "sites" of the potential, while in case (5), the propa- 
gation constant is sufficiently far away from the band edge, and the solitary wave is well 
localized. (The profiles of the corresponding solutions are similar to those of the top and 
bottom solutions shown in Fig. 3 of [Ij.) Case (c) is that of the nonlinear Schrodinger 
equation in two spatial dimensions. In all cases, we apply three methods: the generalized 
Petviashvili method (j2.10p without any acceleration, the same method with the Aitken's 
acceleration (jS.lip performed after every third iteration (naccei = 3), and the mode elimi- 
nation method (|3.7p . (j3.8p with various values of s (see the paragraph including Eq. ()3.9p ). 
The initial condition in all cases is uq = 1.5 e"^^'^"*"^^'*, and the step size At = 1. 

In case (a), the generalized Petviashvili method (j2.10p takes about 950 iterations to 
converge to the accuracy of 10^^''. When the mode elimination technique is applied, starting 
at the moment when the error becomes less or equal to some small value (we chose 10"'^), the 
convergence occurs in about 180 iterations, i.e. more than five times as fast. The evolution 
of the error is shown in Fig. [5] by the thick solid line for the choice s = 1; for smaller values 
of s up to 0.4 which we tried, the error evolution is similar (and the convergence is slightly 
faster). The characteristic feature of this error evolution is that it is nonmonotonic and 
rather irregular. This irregularity is somewhat abated for s < 1, in agreement with our 
discussion in Section 3. Now, when we attempted to apply the Aitken's acceleration to 
the generalized Petviashvili method, we observed quick divergence of the so "accelerated" 
method. We actually tried various values of naccei and At but were unable to make the 
iterations converge. The reason for this is explained at the end of Section 3. In fact, by 
monitoring the error n„ at every iteration, we observed that it contains many nonlocalized 
modes, so that the condition (j3.10p of applicability of the Aitken's acceleration is clearly 
violated in this case. 
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The corresponding results for case (6) are also shown in Fig. [5j There, the mode elimina- 
tion technique accelerates the convergence of the generalized Petviashvili method by about 
a factor of four. The error evolution is much smoother than in case (a). This appears to be 
correlated with the fact, which follows from our monitoring of the error, that the latter is 
dominated by a single eigenmode. Consequently, condition (|3.1U|) is now satisfied, and the 
Steffensen's method (i.e., the generalized Petviashvili method with Aitken's acceleration) 
also converges; see the dotted line in Fig. [5j Let us note that the irregular behavior of the 
error of the Steffensen's method at low values of the error leads to a rather high sensitivity 
of the total number of iterations to the initial condition. For example, we verified that if the 
acceleration is started when the error reaches 10~^ instead of 10^^, the Steffensen's method 
converges to the accuracy of 10^^*^ in about 30 iterations. 

The error evolutions for case (c) are shown in Fig.O The convergence acceleration in this 
case (as, actually, also in case (6)) is not of practical importance because the convergence of 
the non-accelerated generalized Petviashvili method (j2.10p is quite fast (see the thin solid 
line in Fig. [6|). Therefore, below we discuss the results for this method accelerated by the 
mode elimination technique for the sole purpose of highlighting this technique's dependence 
on the parameter s. The error evolution of method (|3.7p . (13. Sp with s = 1, where the 
acceleration is started when the error becomes less or equal to 10"'^, is very irregular (see 
the thick solid line in Fig. [6]), and as a result, the accelerated method takes more iterations 
to converge than the non-accelerated one. Moreover, the evolution of the error also strongly 
depends on the initial condition and on when the acceleration is started. For example, when 
we began the acceleration at the moment of the error reaching 10^^ or 10~^, rather than 
IQ-^, the convergence occurred in about 190 or 100 iterations, respectively. In both cases, 
the error evolution curves were irregular, with several "ups and downs". However, when 
we used values 0.4 < s < 0.8 instead of s = 1, the behavior of the accelerated iterations 
greatly improved. The optimal case of s = 0.7 is shown in Fig. [6] by the medium solid line. 
Both the sensitivity to the "starting moment" of the acceleration and the irregularity of 
the error evolution are suppressed for s < 1, in agreement with the discussion in Section 
3. We also applied the Steffensen's method to this case and found it to converge in about 
the same number of iterations as the mode elimination method with the optinal s; see the 
dotted line in Fig. [H 

Example 4.3 In this and the following Examples, we show that the mode elimination 
technique can be used to accelerate convergence of other iterative methods. In this Example, 
we apply this technique to the squared-operator method (j4.4p . which can converge [8j to 
any given nonfundamental solitary wave of the underlying stationary wave equation. It 
should be noted that in [8] , the efficiency of the so accelerated squared-operator methods 
(referred to there as modified squared-operator methods) was amply demonstrated for a 
number of single and coupled stationary wave equations, both Hamiltonian and dissipative. 
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In all simulations reported in |8], the value of the parameter s in (j3.8p was taken to equal 1. 
Therefore, below we will focus on the dependence of the error evolution on the parameter 
s. 

We apply the squared-operator methods with and without mode elimination to finding 
the second nonfundamental solution of Eq. (j4.ip . This solution for /i = 3 is shown in Fig. [71 
In all cases considered below, we used the initial condition uq = (1 — 2r^) e~^^, = + 
and the step size At = 0.3 (nearly optimal). As the method without mode elimination, we 
used (j4.4p . The method with mode elimination is then a straightforward modification of 
methods dSZD, dSS]) and g31): 

Un+l -Un= (4.7) 



t-It iST-lT r ^Lou)n) {4>slow,n, {LN ^ Lou)n) 

I AT \ -Islow,?! /, 



{N LN Lou)n - Tn -, — r U„ - Fslow, n 77 TTT ^ (Pslow,n 

where, similarly to (j3.8p : 



At, 



LN L(/>slow,n 

t^slow, n — Un Un—1, -L slow, n — ^ , a ' ^slow, n 



•^slow, tiAt (0s1ow, ni -^'/'slow, n) 

(4.8) 

In both cases, with and without mode elimination, we found empirically that the methods 
with r„ given by (14.51) require the initial condition to be closer to the exact solution than 
do the corresponding methods with r„ = 0. On the other hand, the former methods were 
significantly faster than those with r„ = 0. Therefore, we initially used methods (j4.4p or 
(j4.7p with r„ = 0, and when the error reached a small value (we chose 5 • 10~^), switched r„ 
to the expression (j4.5p . The corresponding error evolutions for the accelerated method ()4.7p 
with s = 1 and s = 0.7 (optimal) are shown by the thick and medium lines, while for the 
non-accelerated method (|4.4p without mode elimination, the error evolution is shown by the 
thin line. Note that the behavior of the accelerated method with s < 1 compared to that 
behavior with s = 1 follows the same trends as observed in Example 4.2. Namely, the error 
evolution for the schemes with s < 1 is smoother and much less sensitive to the moment 
when the acceleration starts. Overall, the mode elimination is found to accelerate the 
convergence by a factor between three and four, depending on the choice of the parameter 
s. Finally, we note that the Steffensen's method in this case does not converge. 

Example 4.4 In this last Example, we show that the convergence acceleration technique 
based on mode elimination can also be applied to the ITEM. Here we chose to present 
the results for the version of this method (j2.15p - (j2.17p with power normalization, but the 
technique can be used as well for the ITEM with amplitude normalization [6j. 

For the stationary wave equation (|2.ip written in an equivalent form: 

Lqu = Lqqu — fiu = 0, (4.9) 
the ITEM (I2.15p - (|2.17p with mode elimination can be written as follows: 



Un+l 



p 



{Un+l,Un+l) 



Un+l, (4.10) 
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Un+1 



{Lou) 

^ [^OU)n - Tslow, n TT ^^7 T <Pslow, n 

\'f'slow,n) ^'Pslow, n/ 



At. (4.11) 



Here K is a positive definite self-adjoint operator with constant coefficients (as, e.g., in 

(T \ T {LooUn,K ^Un) , . 

and 

/ ^ S ('/'slow, m -^''/'slow, n) 

'?'slow,n 'I'n— 1) Tslow, n ^ ' T ; '3^slow,n TT TT7 T • 

(^s\ovi,n^'^ \'?'slow,n) 'Pslow, n/ 

(4.13) 

We apply the methods without and with mode elimination — (|2.15p - (|2.17p and (j4.10p - 
(|4.13p . respectively, — to Eq. (j4.6p with Vq = 4 and P = 1.94, whose solution looks similar 
to the top solution in Fig. 3 of [l]. The corresponding propagation constant ^ = 5.01 is 
close to the bandgap edge, and the ITEM without mode elimination converges slowly; see 
the thin line in Fig. O In all simulations, we took Ar = 1 and the operator K of the 
form (j2.12p with c = 1, which yielded the (nearly) optimal convergence rate of the ITEM 
()2.15p - ()2.17p . The error evolutions for the ITEM (j4.10p - (|4.13p with mode elimination are 
shown in Fig. [9] by the thick and medium lines. As in Examples 4.2 and 4.3, the scheme 
with mode elimination provides a severalfold improvement to the convergence rate of the 
ITEM. Also as in those Examples, the error evolution with s < 1 is more regular than that 
with s = 1. 

Thus, from the last three examples, we conclude that in those cases when the iterations 
converge slowly and their acceleration is highly desirable, the mode elimination method 
provides a considerable improvement of the convergence rate (by a factor of several times) . 
Taking s < 1, so that only part of the mode (ti„ — Un-i) would be eliminated, usually results 
in smoother convergence; however, the choice s = 1 still yields a considerable improvement 
of the convergence rate in comparison with that of the non-accelerated iteration method. For 
these slowly convergent cases, the Steff'ensen's method, based on the Aitken's acceleration, 
often diverges. 

Remark In those cases when the step size At is nearly optimal, the error is expected 
to be dominated by two eigenmodes, corresponding to Amax and Amin, since 

(1 + A^axAr) « -(1 + A„,i„AT) (4.14) 

for this At (see (j2.4p and (j2.6p ). Then it seems logical that one would need to eliminate 
both of these eigenmodes, which are proportional to {un — Un-2) and {un — 2un-i + Un-2), 
respectively. We found, however, that although this does result in a smoother error evolution 
than the elimination of just the single mode (n„ — Un-i), it does not yield any consistent 
improvement of the convergence rate compared to the latter case. 
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5 Summary 



In this work, we obtained the following results. 

In Section 2, we compared the linearized operators of the generalized Petviashvili method 
and the ITEM with power normalization. In particular, we showed that while the "primor- 
dial" part of those operators has the same form (|2.19|) . their nonlocal parts involving the 
inner products are different, leading to the eigenvalues of the corresponding operators being 
different. In our simulations we observed that the generalized Petviashvili method converges 
to fundamental solitary waves faster than does the ITEM (although we could prove this 
rigorously only for equations with power-law nonlinearity) . On the other hand, in those 
(rare) cases when both methods converge to a nonfundamental solitary wave, we produced 
an explicit example where the ITEM is faster. 

In Section 3, we proposed a new technique, which we referred to as the mode elimina- 
tion. One application of this technique is that it can obtain nonfundamental solitary waves, 
for which the generalized Petviashvili method would otherwise diverge. The corresponding 
iteration scheme is given by Eqs. (j3.3p and ()3.4p . In Example 4.1 in Section 4, we demon- 
strated that this technique can be superior to an alternative, squared-operator, technique 
[8] when applied to finding lowest- order nonfundamental solutions. However, for finding 
higher-order solutions, the technique of Ref. [8] appears to be more practical. 

As a more important application for the mode elimination technique, we showed that 
it can accelerate the convergence of various iteration methods. This acceleration is most 
significant (by a factor of several times) in those cases when it is most needed, i.e., when 
the convergence of the non-accelerated method is slow. The iteration schemes implementing 
this technique are: Eqs. ()3.7p . ()3.8p for the generalized Petviashvili method; Eqs. ()4.7p . 
(148]) for a squared-operator method (see also Ref. [8]); and Eqs. (f4T0]l - (l4l3]l for the ITEM 
with power normalization. 
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Figure 1: Schematics of the spectra of the operators found on the r.h.s.'es of ()2.19p (a), 
(j2.14p (b), and (j2.18p (c). The circles show the location of discrete eigenvalues. The cross 
on the right of panel (b) indicates the disappearance of the eigenvalue compared to panel 
(a) . The thick dashed line in panel (c) shows a sample function of Eq. (j2.2ip . The left edge 
of the continuous spectrum is located at A = —1. It is assumed that {u, ip^) = 0; see text 
after I^Mi)- 
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Figure 2: (a): The anti-symmetric solution of Eq. ()2.25p with fj, = 1.43 {P = 10). 
(b): Error evolutions, starting with an anti-symmetric initial condition, for the generahzed 
Petviashvih method (solid line) the optimally accelerated ITEM (dashed line). 
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Figure 3: The actual spectra of the linearized operator L of Eq. (j2.25p (a) and of the 
operators N^'^L, where N is given by (j2.12p . with c = 1.5 (b) and c = 5.04 (c). The 
operator in (b) has a very short interval of continuous spectrum between —1 and —0.95. 
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Figure 4: The first nonfundamental solution of Eq. ()4.ip witli /x = 8. 
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Figure 5: The evolution of the error in cases (a) and (b) of Example 4.2. Thin solid: non- 
accelerated method (I2.10p for case (a); thick solid: method ()3.7p . ()3.8p with s = 1 for case 
(a); thin dashed: non-accelerated method (j2.10p for case (b); thick dashed: method (j3.7p . 
(jS.Sp with s = 1 for case (b); thick dotted line: StefFensen's method for case (b). 
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Figure 6: The evolution of the error in case (c) of Example 4.2. Thin solid: non-accelerated 
method (|2.10p : thick sohd: method (|3.7p . (|3.8|) with s = 1; medium solid: method (|3.7p . 
(|3.8p with s = 0.7; thick dotted line: Steffensen's method. 




Figure 7: The second nonfundamental solution of Eq. (|4.1|) with /u = 3. 
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Figure 8: The evolution of the error in Example 4.3; in all cases, the application of ()4.5p 
and/or acceleration is begun when the error norm reaches 5 • 10~^. Thin line: method ()4.4p 
(no mode elimination); thick line: method (|4.7|) with s = 1; medium line: method (|4.7p 
with s = 0.7. 
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Figure 9: The evolution of the error in Example 4.4. Thin line: optimally accelerated 
(with respect to parameter c in operator K, see text) ITEM (|2.15p ~ ()2.17p without mode 
elimination. Thick and medium lines: ITEM (j4.10p - (j4.13p with mode elimination with 
s = 1 (thick) and s = 0.7 (medium). The application of mode elimination begins at the 
first iteration. 
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